‘Easy’ exercises

4E1

The likelihood is

\[ y_i \sim \mathcal{N}(\mu, \sigma) \]

4E2

Two parameters are in the posterior distribution: \(\mu\) and \(\sigma\).

4E3

\[ \mathrm{P}(\mu, \sigma \mid y) = \frac{\prod_i\mathcal{N}(y_i \mid \mu, \sigma )\mathcal{N}(\mu \mid 0, 10)\mathrm{Exp}(\sigma \mid 1)}{\int\int\prod_i\mathcal{N}(y_i \mid \mu, \sigma )\mathcal{N}(\mu \mid 0, 10)\mathrm{Exp}(\sigma \mid 1) d\mu d\sigma} \]

4E4

The linear model is \(\mu_i = \alpha + \beta x_i\).

4E5

Three parameters: \(\alpha\), \(\beta\) and \(\sigma\). \(\mu\) is now determined by \(\alpha\) and \(\beta\).

‘Medium’ Exercises

4M1

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.1     ✓ dplyr   1.0.2
## ✓ tidyr   1.1.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

4M2

4M3

\[ y \sim \mathcal{N}(\mu, \sigma) \\ \mu = \alpha + \beta x \\ \alpha \sim \mathcal{N}(0, 10) \\ \beta \sim \mathcal{U}(0, 1) \\ \sigma \sim \mathrm{Exp}(1) \]

4M4

We choose \(\alpha\) to be a reasonable distribution for adults, and \(\beta\) to be a reasonable distibution of slope centered on zero, because we have no reason to believe right now that height is related to year.

\[ h_i \sim \mathcal{N}(\mu_i, \sigma) \\ \mu_i = \alpha + \beta(x_i - \bar{x}) \\ \alpha \sim \mathcal{N}(170, 20) \\ \beta \sim \mathcal{N}(0, 1) \\ \sigma \sim \mathrm{Exp}(1) \]

4M5

If every student got taller each year, then the mean height would get taller each year, which means that \(\beta\) would have some positive slope. So we would adjust \(\beta\) to be centered on a number greater than zero.

4M6

If the variance is less than 64 then the standard deviation is less than 8. Let’s take a look at the current distribution for \(sigma\).

OK, so this seems to discount any \(\sigma\) > 4. We might want to adjust our rate a bit:

That seems a bit better, so would probably change my rate in the exponential distribution to \(0.5\).

4M7

## Loading required package: rstan
## Loading required package: StanHeaders
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:tidyr':
## 
##     extract
## Loading required package: parallel
## rethinking (Version 2.13)
## 
## Attaching package: 'rethinking'
## The following object is masked from 'package:purrr':
## 
##     map
## The following object is masked from 'package:stats':
## 
##     rstudent
##                   a             b         sigma
## a      7.306605e-02 -3.833504e-08  6.164588e-05
## b     -3.833504e-08  1.754883e-03 -2.272698e-05
## sigma  6.164588e-05 -2.272698e-05  3.653990e-02
##                  a             b         sigma
## a      3.596110602 -0.0783192821  0.0092321087
## b     -0.078319282  0.0017410942 -0.0002015149
## sigma  0.009232109 -0.0002015149  0.0365751797

Centering around the mean has almost eliminated the covariances. Let’s compare posterior models:

The posteriors are pretty much identical whether weight-centered or not.

4M8

It appears that the combination of knots and weight width influences the general wiggliness of the fit lines.

‘Hard’ exercises

4H1

##   weight predicted_height      lwr      upr
## 1  46.95         156.0764 148.2424 164.4932
## 2  43.72         153.6290 145.7414 161.2421
## 3  64.78         178.8947 170.4528 186.9866
## 4  32.59         143.3472 135.4036 150.8580
## 5  54.63         163.0123 155.5557 170.7851

4H2

Let’s get the data first:

Now let’s define a linear model - we choose our priors on the basis that we have no idea about a child’s height:

Let’s get an estimate for b:

##             mean         sd       5.5%      94.5%
## a     189.167000 2.13872602 185.748902 192.585097
## b       2.686824 0.06811439   2.577964   2.795684
## sigma   8.446547 0.43189714   7.756292   9.136802

For every 10 units of weight, the model predicts an additional height of around 27cm. Let’s chart this model.

Weights seems to have a linear relationship with height in the mid-range of approx 15-35 weight units, but the model underestimates height for this range. This is because the model overestimates height at the extreme ends. It looks like a parabolic model might be a better fit here. I would try a quadratic term next.

4H3

This is a remarkably good fit. So each multiple of weight by \(e\) results in an increase in height of approx 49cm. This would explain the shape of the fit. Below we estimate the model and plot it.

4H4

Let’s simulate our prior for various values of our parameters:

Many of these priors are way outside the boundaries of reality. Now we learned from the previous question that when the relationship is linear on the log, the intercept is negative let’s try a negative value for \(\alpha\):

Hasn’t helped a lot, the spread is very wide, so we will reduce the spread of \(\alpha\) and of the other coefficients:

This is better. But we need to make sure that the curves all move upwards. So we should change the distribution for \(\beta_1\) to a plain normal distribution to restrict it more. After playing around with a lot of the coefficient distributions we get to this, which seems a reasonable prior:

4H5

First we try a simple linear model:

Let’s simulate the MAP, 97% confidence interval for \(\mu\) and 97% prediction interval.

And graph the results:

This suggests that higher temperature influences earlier blooming. Now we try a quadratic model:

This suggests no major relationship. Let’s try a cubic:

This suggests that there is a range that does not affect first bloosom day, but that highjer temperatures outside this range relate to earlier blossoming.

Finally lets try a splines model based on cubic basis functions:

Let’s look at a similar relationship between year and March temperature:

## Caution, model may not have converged.
## Code 1: Maximum iterations reached.

Finally let’s pull a previous analysis so we can line these up nicely in a way that suggests possible temperature causality: